Specifically, we will explore the results of 438 datapoints from 15 studies conducted in 3 oceans on 36 species within 23 genera.
## n
## 1 416
## Ref n
## 1 DS34 157
## 2 DS04 58
## 3 DS68 48
## 4 DS16 46
## 5 DS17 24
## 6 DS40 20
## 7 DS10 16
## 8 DS03 12
## 9 DS43 12
## 10 DS20 8
## 11 DS12 7
## 12 DS24 4
## 13 DS28 4
## Ref_name n
## 1 Weber et al. (2006) 157
## 2 Duckworth et al. (2017) 58
## 3 Flores et al. (2012) 48
## 4 Philipp and Fabricius (2003) 46
## 5 Piniak (2007) 24
## 6 Abdel-Salam (1989) PhD Chapter 3 20
## 7 Junjie et al. (2014) 16
## 8 Bessell-Browne et al. (2017a) 12
## 9 Peters and Pilson (1985) 12
## 10 Riegl and Branch (1995) 8
## 11 Loiola et al. (2013) 7
## 12 Sheridan et al. (2014) 4
## 13 Sofonia and Anthony (2008) 4
## Ocean n
## 1 Pacific 349
## 2 Atlantic 39
## 3 Indo-Pacific 16
## 4 Indian 12
## Selecting by n
## Gsp n
## 1 Acropora_millepora 28
## 2 Astrangia_poculata 12
## 3 Galaxea_fascicularis 11
## 4 Goniopora_somaliensis 8
## 5 Montipora_aequituberculata 24
## 6 Montipora_capitata 12
## 7 Montipora_capricornis 20
## 8 Montipora_peltiformis 167
## 9 Porites_lobata/lutea 39
## 10 Turbinaria_reniformis 25
## Selecting by n
## Updated_Genus n
## 1 Montipora 236
## 2 Porites 43
## 3 Turbinaria 32
## 4 Acropora 30
## 5 Astrangia 12
## 6 Galaxea 11
## 7 Goniopora 8
## 8 Mussismilia 7
## 9 Orbicella 6
## 10 Ctenactis 3
## 11 Echinopora 3
## 12 Merulina 3
## 13 Pachyseris 3
## 14 Pectinia 3
Specifically, we will explore the results of 251 datapoints from 7 studies conducted in 3 oceans on 9 species within 6 genera (shown below only the top 10 species and genera, by number of datapoints).
## n
## 1 251
## Ref n
## 1 SS06 146
## 2 SS08 46
## 3 SS17b 20
## 4 SS17a 17
## 5 SS05 12
## 6 SS04 9
## 7 SS25 1
## Ref_name n
## 1 Browne et al. 2015 146
## 2 Flores et al. 2012 46
## 3 Liu et al. 2015 37
## 4 Browne et al. 2014 12
## 5 Bessell-Browne et al. 2017 9
## 6 Te 2001 Chapter 6 1
## Ocean n
## 1 Indian 146
## 2 Pacific 93
## 3 Indo-Pacific (e.g. Singapore) 12
## Gsp n
## 1 Platygyra_sinensis 54
## 2 Merulina_ampliata 53
## 3 Pachyseris_speciosa 51
## 4 Acropora_muricata 37
## 5 Montipora_aequituberculata 30
## 6 Acropora_millepora 19
## 7 Montipora_capricornis 3
## 8 Porites_lutea, lobata 3
## 9 Montipora_verrucosa 1
## Updated_Genus n
## 1 Acropora 56
## 2 Platygyra 54
## 3 Merulina 53
## 4 Pachyseris 51
## 5 Montipora 34
## 6 Porites 3
First let’s explore all data from all species for which there is binary data about the presence of a reduction in ‘photosynthetic efficiency’ or ‘P/R ratio’ as a result of exposure to deposited sediment.
Now let’s calculate the thresholds for each category of photosynthesis, based on the binary data explored above.
## LOAEL
## 1 25
## NOAEL
## 1 25
## LOAEL
## 1 0.5
## NOAEL
## 1 0.5
## LOAEL
## 1 26.4
## NOAEL
## 1 26.4
## LOAEL
## 1 2
## NOAEL
## 1 2
## LOAEL
## 1 35.8
## NOAEL
## 1 35.8
## LOAEL
## 1 56
## NOAEL
## 1 56
## LOAEL
## 1 35.8
## NOAEL
## 1 35.8
## LOAEL
## 1 0.083
## NOAEL
## 1 0.083
## n
## 1 249
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## Data: photo_effDS2
## Models:
## modDS6: Binary_reduced_FvFm ~ Sed_level_stand_mg * Sed_exposure_days +
## modDS6: (1 | Gsp)
## modDS9: Binary_reduced_FvFm ~ Sed_level_stand_mg * Sed_exposure_days +
## modDS9: (1 | Ref)
## modDS12: Binary_reduced_FvFm ~ Sed_level_stand_mg * Sed_exposure_days +
## modDS12: (1 | Gsp) + (1 | Ref)
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## modDS6 5 327.12 344.71 -158.56 317.12
## modDS9 5 314.54 332.13 -152.27 304.54 12.5790 0 <2e-16 ***
## modDS12 6 314.16 335.26 -151.08 302.16 2.3866 1 0.1224
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## Data: photo_effDS2
## Models:
## modDS7: Binary_reduced_FvFm ~ Sed_level_stand_mg + (1 | Ref)
## modDS8: Binary_reduced_FvFm ~ Sed_level_stand_mg + Sed_exposure_days +
## modDS8: (1 | Ref)
## modDS9: Binary_reduced_FvFm ~ Sed_level_stand_mg * Sed_exposure_days +
## modDS9: (1 | Ref)
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## modDS7 3 312.94 323.49 -153.47 306.94
## modDS8 4 313.04 327.11 -152.52 305.04 1.8973 1 0.1684
## modDS9 5 314.54 332.13 -152.27 304.54 0.4953 1 0.4816
## Data: photo_effDS2
## Models:
## modDS10: Binary_reduced_FvFm ~ Sed_level_stand_mg + (1 | Gsp) + (1 | Ref)
## modDS11: Binary_reduced_FvFm ~ Sed_level_stand_mg + Sed_exposure_days +
## modDS11: (1 | Gsp) + (1 | Ref)
## modDS12: Binary_reduced_FvFm ~ Sed_level_stand_mg * Sed_exposure_days +
## modDS12: (1 | Gsp) + (1 | Ref)
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## modDS10 4 313.93 328.00 -152.97 305.93
## modDS11 5 312.74 330.33 -151.37 302.74 3.1894 1 0.07412 .
## modDS12 6 314.16 335.26 -151.08 302.16 0.5867 1 0.44371
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## Data: photo_effDS2
## Models:
## modDS4: Binary_reduced_FvFm ~ Sed_level_stand_mg + (1 | Gsp)
## modDS7: Binary_reduced_FvFm ~ Sed_level_stand_mg + (1 | Ref)
## modDS10: Binary_reduced_FvFm ~ Sed_level_stand_mg + (1 | Gsp) + (1 | Ref)
## modDS11: Binary_reduced_FvFm ~ Sed_level_stand_mg + Sed_exposure_days +
## modDS11: (1 | Gsp) + (1 | Ref)
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## modDS4 3 328.68 339.23 -161.34 322.68
## modDS7 3 312.94 323.49 -153.47 306.94 15.7408 0 < 2e-16 ***
## modDS10 4 313.93 328.00 -152.97 305.93 1.0031 1 0.31656
## modDS11 5 312.74 330.33 -151.37 302.74 3.1894 1 0.07412 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: binomial ( logit )
## Formula: Binary_reduced_FvFm ~ Sed_level_stand_mg + (1 | Ref)
## Data: photo_effDS2
##
## AIC BIC logLik deviance df.resid
## 312.9 323.5 -153.5 306.9 246
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.8252 -0.8127 -0.3036 1.0864 3.3269
##
## Random effects:
## Groups Name Variance Std.Dev.
## Ref (Intercept) 2.737 1.655
## Number of obs: 249, groups: Ref, 9
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -1.0456482 0.6970200 -1.500 0.134
## Sed_level_stand_mg -0.0008905 0.0022947 -0.388 0.698
## [1] 0.6706827
##
## p 0 1
## 0 148 76
## 1 6 19
## Setting levels: control = 0, case = 1
## Setting direction: controls < cases
## Area under the curve: 0.7535
## R2m R2c
## theoretical 0.001289846 0.4548780
## delta 0.001055897 0.3723734
## $Ref
## Groups Name Std.Dev.
## Ref (Intercept) 1.6545
## [1] 5.230593
There is no significant relationship between deposited sediment exposure and the odds of reduced photosynthetic efficiency (GLMM z = -0.388, p = 0.698).
#Overlaying glmm results on figure, but prediction line is jagged...
pred_modDS10 <- predict(modDS10, newdata=photo_effDS2, type="response")
photo_effDS3 <- cbind(photo_effDS2, pred_modDS10)
photo_effDS_plot <- ggplot(data = photo_effDS3) +
geom_point(mapping = aes(
x = Sed_level_stand_mg,
y = Binary_reduced_FvFm)) +
labs(x = expression("Sediment exposure concentration (mg/cm"^"2"*"/day)"),
y = "Reduced photosynthetic efficiency\ndue to sediment exposure?") +
scale_x_log10(limits=c(0.01,1000),breaks=c(0.01,0.1,1,10,100,1000),
label=c("0.01","0.1","1","10","100","1000")) +
geom_line(aes(x = Sed_level_stand_mg, y = pred_modDS10), inherit.aes=FALSE) +
theme_bw()
photo_effDS_plot
That plot is confusing to interpret. Let’s see if I can plot the average marginal probability, i.e., the average change in probability of the outcome across the range of the predictor of interest. This is described in some detail at the following, useful website: https://stats.idre.ucla.edu/r/dae/mixed-effects-logistic-regression/
jvalues <- with(photo_effDS2, seq(from = min(Sed_level_stand_mg), to = max(Sed_level_stand_mg), length.out = 100))
# calculate predicted probabilities and store in a list
pp <- lapply(jvalues, function(j) {
photo_effDS2$Sed_level_stand_mg <- j
predict(modDS10, newdata = photo_effDS2, type = "response")
})
# get the means with lower and upper quartiles
plotdat <- t(sapply(pp, function(x) {
c(M = mean(x), med = median(x), quantile(x, c(0.25, 0.75)))
}))
# add in Sed_level_stand_mg values and convert to data frame
plotdat <- as.data.frame(cbind(plotdat, jvalues))
# better names and show the first few rows
colnames(plotdat) <- c("PredictedProbability", "MedianProbability",
"LowerQuantile", "UpperQuantile", "Sed_conc")
#head(plotdat)
# plot average marginal predicted probabilities
ggplot(plotdat, aes(x = Sed_conc)) +
geom_line(aes(y = PredictedProbability), size = 2) +
geom_line(aes(y = MedianProbability), size = 0.5) +
ylim(c(0, 1))
ggplot(plotdat, aes(x = Sed_conc)) +
geom_ribbon(aes(ymin = LowerQuantile, ymax = UpperQuantile), alpha = 0.15) +
geom_line(aes(y = PredictedProbability), size = 2) +
geom_line(aes(y = MedianProbability), size = 0.5) +
ylim(c(0, 1))
#Overlaying average marginal predicted probabilities on figure
#by Ref
photo_effDS_plot2 <- ggplot() +
geom_point(data = photo_effDS2, mapping = aes(
x = Sed_level_stand_mg,
y = Binary_reduced_FvFm,
color=Ref)) +
labs(x = expression("Sediment exposure concentration (mg/cm"^"2"*"/day)"),
y = "Predicted probability of reduced photosynthetic\nefficiency due to sediment exposure",
color = "Binary data\nby Study",
linetype = "Predicted\nProbability") +
scale_x_log10(limits=c(0.01,1000), breaks=c(0.01,0.1,1,10,100,1000,loaelDS1),
label=c("0.01","0.1","1","10","100","1000",round(loaelDS1,digits=1))) +
geom_ribbon(data = plotdat, aes(x = Sed_conc, y = PredictedProbability,
ymin = LowerQuantile, ymax = UpperQuantile),
alpha = 0.15) +
geom_line(data = plotdat, aes(x = Sed_conc, y = PredictedProbability,
linetype = "twodash")) +
geom_line(data = plotdat, aes(x = Sed_conc, y = MedianProbability,
linetype = "solid")) +
geom_vline(xintercept=loaelDS1, linetype="dashed", color = "red") +
theme_bw() +
scale_linetype_manual(values=c("twodash", "solid"), labels = c("Median","Mean"))
photo_effDS_plot2
## n
## 1 25
## Ref n
## 1 DS10 8
## 2 DS24 2
## 3 DS40 9
## 4 DS43 6
## Ref Gsp
## 1 DS10 Galaxea_fascicularis
## 5 DS10 Goniopora_somaliensis
## 9 DS24 Montipora_patula
## 11 DS40 Agaricia_agaricites
## 12 DS40 Porites_porites
## 13 DS40 Porites_astreoides
## 14 DS40 Diploria_strigosa
## 15 DS40 Isophyllia_rigida
## 16 DS40 Orbicella_annularis
## 17 DS40 Acropora_cervicornis
## 18 DS40 Meandrina_meandrites
## 20 DS43 Astrangia_poculata
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## Data: PR_ratioDS2
## Models:
## modDSb6: Binary_reduced_PR_ratio ~ Sed_level_stand_mg * Sed_exposure_days +
## modDSb6: (1 | Gsp)
## modDSb9: Binary_reduced_PR_ratio ~ Sed_level_stand_mg * Sed_exposure_days +
## modDSb9: (1 | Ref)
## modDSb12: Binary_reduced_PR_ratio ~ Sed_level_stand_mg * Sed_exposure_days +
## modDSb12: (1 | Gsp) + (1 | Ref)
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## modDSb6 5 16.530 22.625 -3.2652 6.5304
## modDSb9 5 13.379 19.474 -1.6895 3.3791 3.1513 0 <2e-16 ***
## modDSb12 6 15.379 22.693 -1.6897 3.3793 0.0000 1 1
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## Data: PR_ratioDS2
## Models:
## modDSb7: Binary_reduced_PR_ratio ~ Sed_level_stand_mg + (1 | Ref)
## modDSb8: Binary_reduced_PR_ratio ~ Sed_level_stand_mg + Sed_exposure_days +
## modDSb8: (1 | Ref)
## modDSb9: Binary_reduced_PR_ratio ~ Sed_level_stand_mg * Sed_exposure_days +
## modDSb9: (1 | Ref)
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## modDSb7 3 25.440 29.096 -9.7198 19.4396
## modDSb8 4 12.511 17.386 -2.2554 4.5108 14.9288 1 0.0001116 ***
## modDSb9 5 13.379 19.474 -1.6895 3.3791 1.1317 1 0.2874148
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## Data: PR_ratioDS2
## Models:
## modDSb4: Binary_reduced_PR_ratio ~ Sed_level_stand_mg + (1 | Gsp)
## modDSb7: Binary_reduced_PR_ratio ~ Sed_level_stand_mg + (1 | Ref)
## modDSb10: Binary_reduced_PR_ratio ~ Sed_level_stand_mg + (1 | Gsp) + (1 |
## modDSb10: Ref)
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## modDSb4 3 32.834 36.491 -13.4171 26.834
## modDSb7 3 25.440 29.096 -9.7198 19.440 7.3946 0 <2e-16 ***
## modDSb10 4 27.440 32.315 -9.7198 19.440 0.0000 1 0.9994
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: binomial ( logit )
## Formula: Binary_reduced_PR_ratio ~ Sed_level_stand_mg + Sed_exposure_days +
## (1 | Ref)
## Data: PR_ratioDS2
##
## AIC BIC logLik deviance df.resid
## 12.5 17.4 -2.3 4.5 21
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -0.0005026 -0.0000769 0.0000000 0.0111328 0.0111328
##
## Random effects:
## Groups Name Variance Std.Dev.
## Ref (Intercept) 78355 279.9
## Number of obs: 25, groups: Ref, 4
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -76.1954 20.5559 -3.707 0.000210 ***
## Sed_level_stand_mg -0.0217 0.1567 -0.138 0.889872
## Sed_exposure_days 4.3993 1.1938 3.685 0.000228 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## [1] 1
##
## p 0 1
## 0 12 0
## 1 0 13
## Setting levels: control = 0, case = 1
## Setting direction: controls < cases
## Area under the curve: 1
## R2m R2c
## theoretical 0.01728004 0.9999587
## delta 0.01725835 0.9987037
## $Ref
## Groups Name Std.Dev.
## Ref (Intercept) 279.92
## [1] 3.694971e+121
## Est LL UL
## (Intercept) 1.155860e-23 8.600357e-36 1.553438e-11
## Sed_level_stand_mg 9.850746e-01 7.962031e-01 1.218749e+00
## Sed_exposure_days 2.110169e+01 4.168496e+00 1.068206e+02
For every doubling in exposure duration of deposited sediment, the odds of reduced P/R ratio increases by 21.1 times (95% CI 4.2, 106.8; GLMM z = 3.685, p = 0.0002), after accounting for exposure concentration. There is no significant relationship between sediment exposure concentration and the odds of reduced P/R ratio (GLMM z = -0.138, p = 0.890). However, all the variability in the model that describes this relationship is attributable to study (R2 = 0.017 with fixed effects only, R2 = 1.000 with fixed effects + random effect of study).
#Overlaying glmm results on figure, but prediction line is jagged...
pred_modDSb8 <- predict(modDSb8, newdata=PR_ratioDS2, type="response")
PR_ratioDS3 <- cbind(PR_ratioDS2, pred_modDSb8)
PR_ratioDS_plot <- ggplot(data = PR_ratioDS3) +
geom_point(mapping = aes(
x = Sed_level_stand_mg,
y = Binary_reduced_PR_ratio,
color = Ref)) +
labs(x = expression("Sediment exposure concentration (mg/cm"^"2"*"/day)"),
y = "Reduced P/R ratio due to sediment exposure?",
color = "Study") +
scale_x_continuous(limits=c(0,200)) +
geom_line(aes(x = Sed_level_stand_mg, y = pred_modDSb8), inherit.aes=FALSE) +
theme_bw()
PR_ratioDS_plot
That plot is confusing to interpret. Let’s see if I can plot the average marginal probability, i.e., the average change in probability of the outcome across the range of the predictor of interest. This is described in some detail at the following, useful website: https://stats.idre.ucla.edu/r/dae/mixed-effects-logistic-regression/
jvalues <- with(PR_ratioDS2, seq(from = min(Sed_level_stand_mg), to = max(Sed_level_stand_mg), length.out = 100))
# calculate predicted probabilities and store in a list
pp <- lapply(jvalues, function(j) {
PR_ratioDS2$Sed_level_stand_mg <- j
predict(modDSb8, newdata = PR_ratioDS2, type = "response")
})
# get the means with lower and upper quartiles
plotdat <- t(sapply(pp, function(x) {
c(M = mean(x), med = median(x), quantile(x, c(0.25, 0.75)))
}))
# add in Sed_level_stand_mg values and convert to data frame
plotdat <- as.data.frame(cbind(plotdat, jvalues))
# better names and show the first few rows
colnames(plotdat) <- c("PredictedProbability", "MedianProbability",
"LowerQuantile", "UpperQuantile", "Sed_conc")
#head(plotdat)
# plot average marginal predicted probabilities
ggplot(plotdat, aes(x = Sed_conc)) +
geom_line(aes(y = PredictedProbability), size = 2) +
geom_line(aes(y = MedianProbability), size = 0.5) +
ylim(c(0, 1))
ggplot(plotdat, aes(x = Sed_conc)) +
geom_ribbon(aes(ymin = LowerQuantile, ymax = UpperQuantile), alpha = 0.15) +
geom_line(aes(y = PredictedProbability), size = 2) +
geom_line(aes(y = MedianProbability), size = 0.5) +
ylim(c(0, 1))
#Overlaying average marginal predicted probabilities on figure
#by Ref
PR_ratioDS_plot2 <- ggplot() +
geom_point(data = PR_ratioDS2, mapping = aes(
x = Sed_level_stand_mg,
y = Binary_reduced_PR_ratio,
color = Ref)) +
labs(x = expression("Sediment exposure concentration (mg/cm"^"2"*"/day)"),
y = "Predicted probability of reduced\nP/R ratio due to sediment exposure",
color = "Binary data\nby Study",
linetype = "Predicted\nProbability") +
scale_x_continuous(limits=c(0,200), breaks=c(0,50,100,150,200,loaelDS2),
labels=c("0","50","100","150","200",round(loaelDS2,digits=1))) +
geom_ribbon(data = plotdat, aes(x = Sed_conc, y = PredictedProbability,
ymin = LowerQuantile, ymax = UpperQuantile),
alpha = 0.15) +
geom_line(data = plotdat, aes(x = Sed_conc, y = PredictedProbability,
linetype = "twodash")) +
geom_line(data = plotdat, aes(x = Sed_conc, y = MedianProbability,
linetype = "solid")) +
geom_vline(xintercept=loaelDS2, linetype="dashed", color = "red") +
theme_bw() +
scale_linetype_manual(values=c("twodash", "solid"), labels = c("Median","Mean"))
PR_ratioDS_plot2
jvalues <- with(PR_ratioDS2, seq(from = min(Sed_exposure_days), to = max(Sed_exposure_days), length.out = 100))
# calculate predicted probabilities and store in a list
pp <- lapply(jvalues, function(j) {
PR_ratioDS2$Sed_exposure_days <- j
predict(modDSb8, newdata = PR_ratioDS2, type = "response")
})
# get the means with lower and upper quartiles
plotdat <- t(sapply(pp, function(x) {
c(M = mean(x), med = median(x), quantile(x, c(0.25, 0.75)))
}))
# add in Sed_exposure_days values and convert to data frame
plotdat <- as.data.frame(cbind(plotdat, jvalues))
# better names and show the first few rows
colnames(plotdat) <- c("PredictedProbability", "MedianProbability",
"LowerQuantile", "UpperQuantile", "Sed_dur")
#head(plotdat)
# plot average marginal predicted probabilities
ggplot(plotdat, aes(x = Sed_dur)) +
geom_line(aes(y = PredictedProbability), size = 2) +
geom_line(aes(y = MedianProbability), size = 0.5) +
ylim(c(0, 1))
ggplot(plotdat, aes(x = Sed_dur)) +
geom_ribbon(aes(ymin = LowerQuantile, ymax = UpperQuantile), alpha = 0.15) +
geom_line(aes(y = PredictedProbability), size = 2) +
geom_line(aes(y = MedianProbability), size = 0.5) +
ylim(c(0, 1))
#Overlaying average marginal predicted probabilities on figure
#by Ref
PR_ratioDS_plot2 <- ggplot() +
geom_point(data = PR_ratioDS2, mapping = aes(
x = Sed_exposure_days,
y = Binary_reduced_PR_ratio,
color = Ref)) +
labs(x = expression("Sediment exposure concentration (mg/cm"^"2"*"/day)"),
y = "Predicted probability of reduced\nP/R ratio due to sediment exposure",
color = "Study",
linetype = "Predicted\nProbability") +
scale_x_continuous(limits=c(0,30)) +
geom_ribbon(data = plotdat, aes(x = Sed_dur, y = PredictedProbability,
ymin = LowerQuantile, ymax = UpperQuantile),
alpha = 0.15) +
geom_line(data = plotdat, aes(x = Sed_dur, y = PredictedProbability,
linetype = "twodash")) +
geom_line(data = plotdat, aes(x = Sed_dur, y = MedianProbability,
linetype = "solid")) +
theme_bw() +
scale_linetype_manual(values=c("twodash", "solid"), labels = c("Median","Mean"))
PR_ratioDS_plot2
## n
## 1 180
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## Data: photo_effSS2
## Models:
## modSS6: Binary_reduced_FvFm ~ Sed_level_stand_mg * Sed_exposure_days +
## modSS6: (1 | Gsp)
## modSS9: Binary_reduced_FvFm ~ Sed_level_stand_mg * Sed_exposure_days +
## modSS9: (1 | Ref)
## modSS12: Binary_reduced_FvFm ~ Sed_level_stand_mg * Sed_exposure_days +
## modSS12: (1 | Gsp) + (1 | Ref)
## modSS15: Binary_reduced_FvFm ~ Sed_level_stand_mg * Sed_exposure_days +
## modSS15: (1 | Ref_name/Ref)
## modSS18: Binary_reduced_FvFm ~ Sed_level_stand_mg * Sed_exposure_days +
## modSS18: (1 | Gsp) + (1 | Ref_name)
## modSS21: Binary_reduced_FvFm ~ Sed_level_stand_mg * Sed_exposure_days +
## modSS21: (1 | Gsp) + (1 | Ref_name/Ref)
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## modSS6 5 29.924 45.889 -9.9622 19.924
## modSS9 5 36.257 52.222 -13.1285 26.257 0.0000 0 1.00000
## modSS12 6 31.924 51.082 -9.9622 19.924 6.3326 1 0.01185 *
## modSS15 6 38.256 57.414 -13.1281 26.256 0.0000 0 1.00000
## modSS18 6 31.924 51.082 -9.9622 19.924 6.3316 0 < 2e-16 ***
## modSS21 7 33.924 56.275 -9.9622 19.924 0.0000 1 1.00000
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## Data: photo_effSS2
## Models:
## modSS4: Binary_reduced_FvFm ~ Sed_level_stand_mg + (1 | Gsp)
## modSS5: Binary_reduced_FvFm ~ Sed_level_stand_mg + Sed_exposure_days +
## modSS5: (1 | Gsp)
## modSS6: Binary_reduced_FvFm ~ Sed_level_stand_mg * Sed_exposure_days +
## modSS6: (1 | Gsp)
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## modSS4 3 33.228 42.807 -13.6138 27.228
## modSS5 4 28.037 40.808 -10.0183 20.037 7.1911 1 0.007327 **
## modSS6 5 29.924 45.889 -9.9622 19.924 0.1121 1 0.737798
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## Data: photo_effSS2
## Models:
## modSS4: Binary_reduced_FvFm ~ Sed_level_stand_mg + (1 | Gsp)
## modSS7: Binary_reduced_FvFm ~ Sed_level_stand_mg + (1 | Ref)
## modSS10: Binary_reduced_FvFm ~ Sed_level_stand_mg + (1 | Gsp) + (1 | Ref)
## modSS13: Binary_reduced_FvFm ~ Sed_level_stand_mg + (1 | Ref_name/Ref)
## modSS16: Binary_reduced_FvFm ~ Sed_level_stand_mg + (1 | Gsp) + (1 | Ref_name)
## modSS19: Binary_reduced_FvFm ~ Sed_level_stand_mg + (1 | Gsp) + (1 | Ref_name/Ref)
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## modSS4 3 33.228 42.807 -13.614 27.228
## modSS7 3 36.149 45.728 -15.075 30.149 0.0000 0 1.0000
## modSS10 4 35.228 47.999 -13.614 27.228 2.9215 1 0.0874 .
## modSS13 4 38.149 50.921 -15.075 30.149 0.0000 0 1.0000
## modSS16 4 35.228 47.999 -13.614 27.228 2.9215 0 <2e-16 ***
## modSS19 5 37.228 53.192 -13.614 27.228 0.0000 1 1.0000
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: binomial ( logit )
## Formula: Binary_reduced_FvFm ~ Sed_level_stand_mg + Sed_exposure_days +
## (1 | Gsp)
## Data: photo_effSS2
##
## AIC BIC logLik deviance df.resid
## 28.0 40.8 -10.0 20.0 176
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.20089 -0.00726 -0.00048 -0.00015 2.64336
##
## Random effects:
## Groups Name Variance Std.Dev.
## Gsp (Intercept) 94.68 9.731
## Number of obs: 180, groups: Gsp, 8
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -21.30065 9.94516 -2.142 0.0322 *
## Sed_level_stand_mg 0.02395 0.02453 0.976 0.3289
## Sed_exposure_days 0.13636 0.07874 1.732 0.0833 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## [1] 0.9777778
##
## p 0 1
## 0 176 3
## 1 1 0
## Setting levels: control = 0, case = 1
## Setting direction: controls < cases
## Area under the curve: 0.9859
## R2m R2c
## theoretical 0.09587058 0.9696399
## delta 0.08433380 0.8529564
## $Gsp
## Groups Name Std.Dev.
## Gsp (Intercept) 9.7305
## [1] 16823.41
## Est LL UL
## (Intercept) 3.871372e-07 5.248454e-13 0.2855607
## Sed_level_stand_mg 1.016741e+00 9.834110e-01 1.0512007
## Sed_exposure_days 1.099130e+00 9.876226e-01 1.2232273
There is no significant relationship between exposure concentration of suspended sediment and the odds of reduced photosynthetic efficiency (GLMM z = 0.976, p = 0.329), though there is suggestive, but non-significant evidence that for every doubling of exposure duration that the odds of reduced photosynthetic efficiency increase by 1.1 times (95% CI 1.0, 1.2; GLMM z = 1.732, p = 0.083).
#Overlaying glmm results on figure, but prediction line is jagged...
pred_modSS5 <- predict(modSS5, newdata=photo_effSS2, type="response")
photo_effSS3 <- cbind(photo_effSS2, pred_modSS5)
photo_effSS_plot <- ggplot(data = photo_effSS3) +
geom_point(mapping = aes(
x = Sed_level_stand_mg,
y = Binary_reduced_FvFm,
color = Ref)) +
labs(x = "Sediment exposure concentration (mg/L)",
y = "Reduced photosynthetic efficiency due to sediment exposure",
color = "Study") +
scale_x_log10(limits=c(1,250),breaks=c(0.01,0.1,1,10,100,1000),
label=c("0.01","0.1","1","10","100","1000")) +
geom_line(aes(x = Sed_level_stand_mg, y = pred_modSS5), inherit.aes=FALSE) +
theme_bw()
photo_effSS_plot
That plot is confusing to interpret. Let’s see if I can plot the average marginal probability, i.e., the average change in probability of the outcome across the range of the predictor of interest. This is described in some detail at the following, useful website: https://stats.idre.ucla.edu/r/dae/mixed-effects-logistic-regression/
jvalues <- with(photo_effSS2, seq(from = min(Sed_level_stand_mg), to = max(Sed_level_stand_mg), length.out = 100))
# calculate predicted probabilities and store in a list
pp <- lapply(jvalues, function(j) {
photo_effSS2$Sed_level_stand_mg <- j
predict(modSS5, newdata = photo_effSS2, type = "response")
})
# get the means with lower and upper quartiles
plotdat <- t(sapply(pp, function(x) {
c(M = mean(x), med = median(x), quantile(x, c(0.25, 0.75)))
}))
# add in Sed_level_stand_mg values and convert to data frame
plotdat <- as.data.frame(cbind(plotdat, jvalues))
# better names and show the first few rows
colnames(plotdat) <- c("PredictedProbability", "MedianProbability",
"LowerQuantile", "UpperQuantile", "Sed_conc")
#head(plotdat)
# plot average marginal predicted probabilities
ggplot(plotdat, aes(x = Sed_conc)) +
geom_line(aes(y = PredictedProbability), size = 2) +
geom_line(aes(y = MedianProbability), size = 0.5) +
ylim(c(0, 1))
ggplot(plotdat, aes(x = Sed_conc)) +
geom_ribbon(aes(ymin = LowerQuantile, ymax = UpperQuantile), alpha = 0.15) +
geom_line(aes(y = PredictedProbability), size = 2) +
geom_line(aes(y = MedianProbability), size = 0.5) +
ylim(c(0, 1))
#Overlaying average marginal predicted probabilities on figure
#by Ref
photo_effSS_plot2 <- ggplot() +
geom_point(data = photo_effSS2, mapping = aes(
x = Sed_level_stand_mg,
y = Binary_reduced_FvFm,
color = Ref)) +
labs(x = "Sediment exposure concentration (mg/L)",
y = "Predicted probability of reduced photosynthetic\nefficiency due to sediment exposure",
color = "Binary Data\nby Study",
linetype = "Predicted\nProbability") +
scale_x_log10(limits=c(1,max(photo_effSS2$Sed_level_stand_mg)),
breaks=c(0.01,0.1,1,10,100,1000,loaelSS1),
label=c("0.01","0.1","1","10","100","1000",round(loaelSS1,digits = 1))) +
geom_ribbon(data = plotdat, aes(x = Sed_conc, y = PredictedProbability,
ymin = LowerQuantile, ymax = UpperQuantile),
alpha = 0.15) +
geom_line(data = plotdat, aes(x = Sed_conc, y = PredictedProbability,
linetype = "twodash")) +
geom_line(data = plotdat, aes(x = Sed_conc, y = MedianProbability,
linetype = "solid")) +
geom_vline(xintercept=loaelSS1, linetype="dashed", color = "red") +
theme_bw() +
scale_linetype_manual(values=c("twodash", "solid"), labels = c("Median","Mean"))
photo_effSS_plot2
jvalues <- with(photo_effSS2, seq(from = min(Sed_exposure_days), to = max(Sed_exposure_days), length.out = 100))
# calculate predicted probabilities and store in a list
pp <- lapply(jvalues, function(j) {
photo_effSS2$Sed_exposure_days <- j
predict(modSS5, newdata = photo_effSS2, type = "response")
})
# average marginal predicted probability across a few different sediment levels
#sapply(pp[c(1, 20, 40, 60, 80, 100)], mean)
## [1] 0.1797186 0.1960244 0.2142042 0.2333993 0.2535694 0.2746604
# get the means with lower and upper quartiles
plotdat <- t(sapply(pp, function(x) {
c(M = mean(x), med = median(x), quantile(x, c(0.25, 0.75)))
}))
# add in Sed_exposure_days values and convert to data frame
plotdat <- as.data.frame(cbind(plotdat, jvalues))
# better names and show the first few rows
colnames(plotdat) <- c("PredictedProbability", "MedianProbability",
"LowerQuantile", "UpperQuantile", "Sed_dur")
#head(plotdat)
# plot average marginal predicted probabilities
ggplot(plotdat, aes(x = Sed_dur)) +
geom_line(aes(y = PredictedProbability), size = 2) +
geom_line(aes(y = MedianProbability), size = 0.5) +
ylim(c(0, 1))
ggplot(plotdat, aes(x = Sed_dur)) +
geom_ribbon(aes(ymin = LowerQuantile, ymax = UpperQuantile), alpha = 0.15) +
geom_line(aes(y = PredictedProbability), size = 2) +
geom_line(aes(y = MedianProbability), size = 0.5) +
ylim(c(0, 1))
#Overlaying average marginal predicted probabilities on figure
#by Ref
photo_effSS_plot3 <- ggplot() +
geom_point(data = photo_effSS2, mapping = aes(
x = Sed_exposure_days,
y = Binary_reduced_FvFm,
color = Ref)) +
labs(x = "Sediment exposure duration (days)",
y = "Predicted probability of reduced photosynthetic\nefficiency due to sediment exposure",
color = "Binary Data\nby Study",
linetype = "Predicted\nProbability") +
scale_x_continuous(limits=c(0,max(photo_effSS2$Sed_exposure_days))) +
geom_ribbon(data = plotdat, aes(x = Sed_dur, y = PredictedProbability,
ymin = LowerQuantile, ymax = UpperQuantile),
alpha = 0.15) +
geom_line(data = plotdat, aes(x = Sed_dur, y = PredictedProbability,
linetype = "twodash")) +
geom_line(data = plotdat, aes(x = Sed_dur, y = MedianProbability,
linetype = "solid")) +
theme_bw() +
scale_linetype_manual(values=c("twodash", "solid"), labels = c("Median","Mean"))
photo_effSS_plot3
## n
## 1 34
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## Data: PR_ratioSS2
## Models:
## modSSb6: Binary_reduced_PR_ratio ~ Sed_level_stand_mg * Sed_exposure_days +
## modSSb6: (1 | Gsp)
## modSSb9: Binary_reduced_PR_ratio ~ Sed_level_stand_mg * Sed_exposure_days +
## modSSb9: (1 | Ref)
## modSSb12: Binary_reduced_PR_ratio ~ Sed_level_stand_mg * Sed_exposure_days +
## modSSb12: (1 | Gsp) + (1 | Ref)
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## modSSb6 5 36.328 43.959 -13.164 26.328
## modSSb9 5 36.378 44.010 -13.189 26.378 0.0000 0 1.000
## modSSb12 6 38.328 47.486 -13.164 26.328 0.0506 1 0.822
## Data: PR_ratioSS2
## Models:
## modSSb7: Binary_reduced_PR_ratio ~ Sed_level_stand_mg + (1 | Ref)
## modSSb8: Binary_reduced_PR_ratio ~ Sed_level_stand_mg + Sed_exposure_days +
## modSSb8: (1 | Ref)
## modSSb9: Binary_reduced_PR_ratio ~ Sed_level_stand_mg * Sed_exposure_days +
## modSSb9: (1 | Ref)
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## modSSb7 3 33.577 38.156 -13.788 27.577
## modSSb8 4 35.425 41.531 -13.713 27.425 0.1514 1 0.6972
## modSSb9 5 36.378 44.010 -13.189 26.378 1.0468 1 0.3062
## Data: PR_ratioSS2
## Models:
## modSSb4: Binary_reduced_PR_ratio ~ Sed_level_stand_mg + (1 | Gsp)
## modSSb7: Binary_reduced_PR_ratio ~ Sed_level_stand_mg + (1 | Ref)
## modSSb10: Binary_reduced_PR_ratio ~ Sed_level_stand_mg + (1 | Gsp) + (1 |
## modSSb10: Ref)
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## modSSb4 3 33.549 38.128 -13.774 27.549
## modSSb7 3 33.577 38.156 -13.788 27.577 0.0000 0 1.0000
## modSSb10 4 35.549 41.654 -13.774 27.549 0.0276 1 0.8682
## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: binomial ( logit )
## Formula: Binary_reduced_PR_ratio ~ Sed_level_stand_mg + (1 | Ref)
## Data: PR_ratioSS2
##
## AIC BIC logLik deviance df.resid
## 33.6 38.2 -13.8 27.6 31
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -0.6581 -0.4039 -0.4039 -0.3360 2.9765
##
## Random effects:
## Groups Name Variance Std.Dev.
## Ref (Intercept) 1.256e-16 1.121e-08
## Number of obs: 34, groups: Ref, 3
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -2.414418 0.914690 -2.640 0.0083 **
## Sed_level_stand_mg 0.006506 0.007018 0.927 0.3539
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## convergence code: 0
## boundary (singular) fit: see ?isSingular
## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: binomial ( logit )
## Formula: Binary_reduced_PR_ratio ~ Sed_level_stand_mg + (1 | Gsp)
## Data: PR_ratioSS2
##
## AIC BIC logLik deviance df.resid
## 33.5 38.1 -13.8 27.5 31
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -0.6967 -0.4137 -0.3646 -0.3060 2.8420
##
## Random effects:
## Groups Name Variance Std.Dev.
## Gsp (Intercept) 0.1163 0.341
## Number of obs: 34, groups: Gsp, 4
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -2.468090 1.002623 -2.462 0.0138 *
## Sed_level_stand_mg 0.006610 0.007112 0.929 0.3527
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## [1] 0.8529412
##
## p 0 1
## 0 29 5
## Setting levels: control = 0, case = 1
## Setting direction: controls < cases
## Area under the curve: 0.6655
## R2m R2c
## theoretical 0.04833292 0.08081445
## delta 0.02088213 0.03491570
## $Gsp
## Groups Name Std.Dev.
## Gsp (Intercept) 0.34096
## [1] 1.4063
There is no significant relationship between exposure concentration of suspended sediment and the odds of a reduced P/R ratio (GLMM z = 0.929, p = 0.353), but the best-fit model to explain this relationship explains only 8% of the variance, indicating that more data is needed.
#Overlaying glmm results on figure, but prediction line is jagged...
pred_modSSb4 <- predict(modSSb4, newdata=PR_ratioSS2, type="response")
PR_ratioSS3 <- cbind(PR_ratioSS2, pred_modSSb4)
PR_ratioSS_plot <- ggplot(data = PR_ratioSS3) +
geom_point(mapping = aes(
x = Sed_level_stand_mg,
y = Binary_reduced_PR_ratio,
color = Ref)) +
labs(x = "Sediment exposure concentration (mg/L)",
y = "Reduced P/R ratio due to sediment exposure",
color = "Study") +
scale_x_log10(limits=c(10,1000),breaks=c(0.01,0.1,1,10,100,1000),
label=c("0.01","0.1","1","10","100","1000")) +
geom_line(aes(x = Sed_level_stand_mg, y = pred_modSSb4), inherit.aes=FALSE) +
theme_bw()
PR_ratioSS_plot
That plot is confusing to interpret. Let’s see if I can plot the average marginal probability, i.e., the average change in probability of the outcome across the range of the predictor of interest. This is described in some detail at the following, useful website: https://stats.idre.ucla.edu/r/dae/mixed-effects-logistic-regression/
jvalues <- with(PR_ratioSS2, seq(from = min(Sed_level_stand_mg), to = max(Sed_level_stand_mg), length.out = 100))
# calculate predicted probabilities and store in a list
pp <- lapply(jvalues, function(j) {
PR_ratioSS2$Sed_level_stand_mg <- j
predict(modSSb4, newdata = PR_ratioSS2, type = "response")
})
# get the means with lower and upper quartiles
plotdat <- t(sapply(pp, function(x) {
c(M = mean(x), med = median(x), quantile(x, c(0.25, 0.75)))
}))
# add in Sed_level_stand_mg values and convert to data frame
plotdat <- as.data.frame(cbind(plotdat, jvalues))
# better names and show the first few rows
colnames(plotdat) <- c("PredictedProbability", "MedianProbability",
"LowerQuantile", "UpperQuantile", "Sed_conc")
#head(plotdat)
# plot average marginal predicted probabilities
ggplot(plotdat, aes(x = Sed_conc)) +
geom_line(aes(y = PredictedProbability), size = 2) +
geom_line(aes(y = MedianProbability), size = 0.5) +
ylim(c(0, 1))
ggplot(plotdat, aes(x = Sed_conc)) +
geom_ribbon(aes(ymin = LowerQuantile, ymax = UpperQuantile), alpha = 0.15) +
geom_line(aes(y = PredictedProbability), size = 2) +
geom_line(aes(y = MedianProbability), size = 0.5) +
ylim(c(0, 1))
#Overlaying average marginal predicted probabilities on figure
#by Ref
PR_ratioSS_plot2 <- ggplot() +
geom_point(data = PR_ratioSS2, mapping = aes(
x = Sed_level_stand_mg,
y = Binary_reduced_PR_ratio,
color = Ref)) +
labs(x = "Sediment exposure concentration (mg/L)",
y = "Predicted probability of reduced\nP/R ratio to sediment exposure",
color = "Binary Data\nby Study",
linetype = "Predicted\nProbability") +
scale_x_log10(limits=c(10,max(PR_ratioSS2$Sed_level_stand_mg)),
breaks=c(0.01,0.1,1,10,100,1000,loaelSS2),
label=c("0.01","0.1","1","10","100","1000",round(loaelSS2,digits = 1))) +
geom_ribbon(data = plotdat, aes(x = Sed_conc, y = PredictedProbability,
ymin = LowerQuantile, ymax = UpperQuantile),
alpha = 0.15) +
geom_line(data = plotdat, aes(x = Sed_conc, y = PredictedProbability,
linetype = "twodash")) +
geom_line(data = plotdat, aes(x = Sed_conc, y = MedianProbability,
linetype = "solid")) +
geom_vline(xintercept=loaelSS2, linetype="dashed", color = "red") +
theme_bw() +
scale_linetype_manual(values=c("twodash", "solid"), labels = c("Median","Mean"))
PR_ratioSS_plot2
jvalues <- with(PR_ratioSS2, seq(from = min(Sed_exposure_days), to = max(Sed_exposure_days), length.out = 100))
# calculate predicted probabilities and store in a list
pp <- lapply(jvalues, function(j) {
PR_ratioSS2$Sed_exposure_days <- j
predict(modSSb4, newdata = PR_ratioSS2, type = "response")
})
# average marginal predicted probability across a few different sediment levels
#sapply(pp[c(1, 20, 40, 60, 80, 100)], mean)
## [1] 0.1797186 0.1960244 0.2142042 0.2333993 0.2535694 0.2746604
# get the means with lower and upper quartiles
plotdat <- t(sapply(pp, function(x) {
c(M = mean(x), med = median(x), quantile(x, c(0.25, 0.75)))
}))
# add in Sed_exposure_days values and convert to data frame
plotdat <- as.data.frame(cbind(plotdat, jvalues))
# better names and show the first few rows
colnames(plotdat) <- c("PredictedProbability", "MedianProbability",
"LowerQuantile", "UpperQuantile", "Sed_dur")
#head(plotdat)
# plot average marginal predicted probabilities
ggplot(plotdat, aes(x = Sed_dur)) +
geom_line(aes(y = PredictedProbability), size = 2) +
geom_line(aes(y = MedianProbability), size = 0.5) +
ylim(c(0, 1))
ggplot(plotdat, aes(x = Sed_dur)) +
geom_ribbon(aes(ymin = LowerQuantile, ymax = UpperQuantile), alpha = 0.15) +
geom_line(aes(y = PredictedProbability), size = 2) +
geom_line(aes(y = MedianProbability), size = 0.5) +
ylim(c(0, 1))
#Overlaying average marginal predicted probabilities on figure
#by Ref
PR_ratioSS_plot3 <- ggplot() +
geom_point(data = PR_ratioSS2, mapping = aes(
x = Sed_exposure_days,
y = Binary_reduced_PR_ratio,
color = Ref)) +
labs(x = "Sediment exposure duration (days)",
y = "Predicted probability of reduced\nP/R ratio to sediment exposure",
color = "Binary Data\nby Study",
linetype = "Predicted\nProbability") +
scale_x_continuous(limits=c(0,max(PR_ratioSS2$Sed_exposure_days))) +
geom_ribbon(data = plotdat, aes(x = Sed_dur, y = PredictedProbability,
ymin = LowerQuantile, ymax = UpperQuantile),
alpha = 0.15) +
geom_line(data = plotdat, aes(x = Sed_dur, y = PredictedProbability,
linetype = "twodash")) +
geom_line(data = plotdat, aes(x = Sed_dur, y = MedianProbability,
linetype = "solid")) +
theme_bw() +
scale_linetype_manual(values=c("twodash", "solid"), labels = c("Median","Mean"))
PR_ratioSS_plot3